# File combines data on energy use and GDP from WIOD 2016 release with exposures scores

library(readxl)
library(openxlsx)
library(readxlsb)
library(ggplot2)
library(concordance)
library(dplyr)
library(RColorBrewer)
library(gridExtra)

# Set working directory
path = "G:/My Drive/Research/AI and Energy/Drafting/ERL/Replication"
setwd(path)


################################################################################
################################################################################
# 1. Import and Prep WIOD Data
################################################################################
################################################################################
# Initialize dataframe to store relevant data
df.WIOD <- data.frame(matrix(nrow=0,ncol=5))
colnames(df.WIOD)=c("Year","Code","Description","Total Energy Use (TJ)","Output ($MM)")

# Get sector codes and description
df.categories <- as.data.frame(read_excel("Data/Input/WIOD/Gross10/NAMEA_GEU5610_USA.xls",sheet="Sector"))[c(1:57),]
df.categories$Code <- gsub("^.{0,1}", "", df.categories$Code) # Remove "v" in front of codes
# Get emissions data
df.emissions <- read.xlsx("Data/Input/WIOD/CO2 emissions.xlsx",sheet="USA")

# Loop over years and append yearly data to data.frame
for(yr in c(2000:2014)){
  print(yr)
  # Get energy data
  energy.yr <- read_excel("Data/Input/WIOD/Gross10/NAMEA_GEU5610_USA.xls",sheet=as.character(yr))[c(1:57),c(15)]$TOTAL
  # Get gdp data
  gdp.yr <- c(read_xlsb(paste0("Data/Input/WIOD/WIOT",as.character(yr),"_Nov16_ROW.xlsb"),range=paste0(as.character(yr),"!CYK2359:CYK2414"),col_names=F)$column.1,sum(read_xlsb(paste0("Data/Input/WIOD/WIOT",as.character(yr),"_Nov16_ROW.xlsb"),range=paste0(as.character(yr),"!CYK2359:CYK2414"),col_names=F)$column.1))
  # Get Emissions data
  emissions.yr <- c(df.emissions[c(1:56),as.character(yr)],sum(df.emissions[c(1:56),as.character(yr)]))
  df.tmp <- data.frame("Year"=yr,"Code"=df.categories[,1],"Description"=df.categories[,2],"Total Energy Use (TJ)"=energy.yr,"Output ($MM)"=gdp.yr,"CO2 Emissions (kt)"=emissions.yr)
  df.WIOD <- rbind(df.WIOD,df.tmp)
}
colnames(df.WIOD)=c("Year","Code","Description","Total Energy Use (TJ)","Output ($MM)","CO2 Emissions (kt)")


# Deflate prices to 2017USD
gdp.def = read.csv("Data/Input/A191RD3A086NBEA.csv")
gdp.def['Year'] = unlist(lapply(gdp.def['DATE'],FUN=function(x){as.numeric(substr(x,1,4))}))
gdp.def = gdp.def[,c(2,3)]
colnames(gdp.def) = c("GDP_Def","Year")
gdp.def['GDP_Def'] = gdp.def['GDP_Def']/100
df.WIOD = merge(df.WIOD,gdp.def,by="Year")
df.WIOD['Output ($MM)'] = df.WIOD['Output ($MM)']/df.WIOD['GDP_Def']

# Change units for analysis
df.WIOD['Output ($BB)'] = df.WIOD['Output ($MM)']*1e-3
df.WIOD['Total Energy Use (PJ)'] = df.WIOD['Total Energy Use (TJ)']*1e-3


# Create useful parameters
df.WIOD['Energy Intensity (PJ/$BB)'] = df.WIOD['Total Energy Use (PJ)']/df.WIOD['Output ($BB)']
df.WIOD['Emissions Intensity (kt/$BB)'] = df.WIOD['CO2 Emissions (kt)']/df.WIOD['Output ($BB)']
df.WIOD['Emissions Intensity (kt/PJ)'] = df.WIOD['CO2 Emissions (kt)']/df.WIOD['Total Energy Use (PJ)']






# Export WIOD merged data
write.csv(df.WIOD,file="Data/Output/Energy Intensity by Industry.csv")





################################################################################
################################################################################
# 2. Merge NAICS exposure weights (4-digit) with WIOD ISIC classifications
################################################################################
################################################################################

#############
# Create concordance array between ISIC codes and NAICS
#############

# Load NAICS 4-digit exposure weight data
df.NAICS.exp <- read.csv(file="Data/Output/NAICS4_Occupation_wagebillweighted_exposure.csv")

# Initialize array to store mapping
# Get 2-digit ISIC codes
df.ISIC_codes <- read_excel("Data/Input/CL_ACTIVITY_ISIC4_1.0.xls")
ISIC_codes_2_withletter <- array(unlist(df.ISIC_codes[df.ISIC_codes$`Annotation 1 - Hierarchical Level`=="Hierarchical level 2","Code Value"]))
ISIC_codes_2 <- substr(ISIC_codes_2_withletter,2,3)
# Get 4-digit NAICS codes
NAICS_codes_4 <- substr(unique(df.NAICS.exp$NAICS),1,4)
# Some NAICS codes are either multiples of 4 digit codes (e.g. ending in A1) or have bad match. For some, aggregate to 3 digit when possible. For others, aggregate to 2 digit. For others, no match.
missmatch3 <- c("325","327","332","333","337","423","424","445","484","522","523","531","532")
missmatch2 <- c("45","51")

# Create array using all matches
concord.arr.all <- array(0,dim=c(length(ISIC_codes_2),length(NAICS_codes_4)),dimnames=list(ISIC_codes_2,NAICS_codes_4))

# Fill in array using all matches
missing.naics = c()
for(naics in NAICS_codes_4){
  naics_lookup = naics
  
  # Replace with 2 or 3 digit naics when needed
  if(substr(naics,1,3)%in%missmatch3){naics_lookup = substr(naics,1,3)}
  if(substr(naics,1,2)%in%missmatch2){naics_lookup = substr(naics,1,2)}
  
  # Retreive isic code from concordance
  isic = concord_naics_isic(naics_lookup,"NAICS2017","ISIC4",dest.digit=2,all = TRUE)
  
  if(is.na(isic[[1]]$match[1])){
    print(paste0("No match found for ",naics))
    missing.naics = c(missing.naics,naics)
  }else{
    for(m in c(1:length(isic[[1]]$match))){
      concord.arr.all[isic[[1]]$match[m],naics] = concord.arr.all[isic[[1]]$match[m],naics] + isic[[1]]$weight[m]
    }
  }
}

# Create array using best match
concord.arr.best <- array(0,dim=c(length(ISIC_codes_2),length(NAICS_codes_4)),dimnames=list(ISIC_codes_2,NAICS_codes_4))

# Fill in array using best matches
missing.naics = c()
for(naics in NAICS_codes_4){
  naics_lookup = naics
  
  # Replace with 2 or 3 digit naics when needed
  if(substr(naics,1,3)%in%missmatch3){naics_lookup = substr(naics,1,3)}
  if(substr(naics,1,2)%in%missmatch2){naics_lookup = substr(naics,1,2)}
  
  # Retreive isic code from concordance
  isic = concord_naics_isic(naics_lookup,"NAICS2017","ISIC4",dest.digit=2,all = FALSE)
  
  if(is.na(isic)){ # If no match at 2-digit, throw warning
    missing.naics = c(missing.naics,naics)
    print(paste0("No match found for ",naics))
  }else{
    concord.arr.best[isic,naics] = 1
  }
}


#############
# Aggregate concordance array along ISIC codes to match WIOD codes
#############


# Data.frame matching ISIC codes in WIOD to all 2-digit ISIC codes
wiod_to_isic <- data.frame("ISIC2" = ISIC_codes_2,"WIOD_code" = NA)
letter_look_up = function(letter){
  nums = substr(ISIC_codes_2_withletter[substr(ISIC_codes_2_withletter,1,1)==letter],2,3)
  return(nums)
}
for(i in unique(df.WIOD$Code)[1:56]){
  if(nchar(i)==1){ # If WIOD code is a number
    nums = letter_look_up(i)
    for(n in nums){
      wiod_to_isic[wiod_to_isic$ISIC2==n,"WIOD_code"] = i
    }
  }else if(grepl("_",i)){# If WIOD code spans multiple codes
    if(i=="R_S"){# If R/S
      i_subs = c("R","S")
      for(is in i_subs){
        nums = letter_look_up(is)
        for(n in nums){
          wiod_to_isic[wiod_to_isic$ISIC2==n,"WIOD_code"] = i
        }
      }
    }else{# If multiple numbers
      nums = as.character(c(as.numeric(substr(i,2,3)):as.numeric(substr(i,5,6))))
      for(n in nums){
        wiod_to_isic[wiod_to_isic$ISIC2==n,"WIOD_code"] = i
      }
    }
  }else if(nchar(i)==3){
    wiod_to_isic[wiod_to_isic$ISIC2==substr(i,2,3),"WIOD_code"] = i
  }
}

# Aggregate along WIOD codes
concord.arr.all.WIOD <- concord.arr.best.WIOD <- array(0,dim=c(length(unique(df.WIOD$Code)[1:56]),length(NAICS_codes_4)),dimnames=list(unique(df.WIOD$Code)[1:56],NAICS_codes_4))
for(i in dimnames(concord.arr.best.WIOD)[[1]]){
  isics = wiod_to_isic[wiod_to_isic$WIOD_code==i,"ISIC2"]
  if(length(isics)==1){ # If only one 2-digit code in WIOD ISIC code, look up code
    concord.arr.all.WIOD[i,] <- concord.arr.all[isics,]
    concord.arr.best.WIOD[i,] <- concord.arr.best[isics,]
  }else{ # If multiple 2-digit code in WIOD ISIC code, aggregate over ISIC codes
    concord.arr.all.WIOD[i,] <- apply(concord.arr.all[isics,],2,sum)
    concord.arr.best.WIOD[i,] <- apply(concord.arr.best[isics,],2,sum)
  }
}




##########
# Merge NAICS data with WIOD
##########
# Get list of exposure metrics in NAICS dataset
expmeasures = colnames(df.NAICS.exp)[substr(colnames(df.NAICS.exp),1,8)=="expshare"]

# Duplicate WIOD data for merging
df.WIOD.exp.all <- df.WIOD.exp.best <- df.WIOD
df.WIOD.exp.all[expmeasures] <- df.WIOD.exp.best[expmeasures] <- 0
for(naic in unique(df.NAICS.exp$NAICS)){
  if(sum(concord.arr.best.WIOD[,substr(naic,1,4)])==0){next} # If NAICS doesn't match, skip
  
  # Fill in for all match
  cells.all = names(which(concord.arr.all.WIOD[,substr(naic,1,4)]>0))
  weights.all = concord.arr.all.WIOD[cells.all,substr(naic,1,4)]
  if(length(cells.all)>1){ # If NAICS matches to multiple isic codes, weight using weights from concord function
    for(j in c(1:length(cells.all))){
      for(i in expmeasures){
        df.WIOD.exp.all[df.WIOD.exp.all$Code==cells.all[j],i] = df.WIOD.exp.all[df.WIOD.exp.all$Code==cells.all[j],i] + df.NAICS.exp[df.NAICS.exp$NAICS==naic,i]*weights.all[j]
      }
    }
  }else if(length(cells.all)==1){
    for(i in expmeasures){
      df.WIOD.exp.all[df.WIOD.exp.all$Code==cells.all,i] = df.WIOD.exp.all[df.WIOD.exp.all$Code==cells.all,i] + df.NAICS.exp[df.NAICS.exp$NAICS==naic,i]
    }
  }
  
  # Fill in for best match
  cells.best = names(which(concord.arr.best.WIOD[,substr(naic,1,4)]>0))
  weights.best = concord.arr.best.WIOD[cells.best,substr(naic,1,4)]
  for(i in expmeasures){
    df.WIOD.exp.best[df.WIOD.exp.best$Code==cells.best,i] = df.WIOD.exp.best[df.WIOD.exp.best$Code==cells.best,i] + df.NAICS.exp[df.NAICS.exp$NAICS==naic,i]
  }
}

# Adjust to be share of VA_k rather than relative to GDP
for(i in expmeasures[c(2,4,6,7,10,12,14,16)]){
  df.WIOD.exp.best[i] = df.WIOD.exp.best[i]*23e9/(df.WIOD.exp.best['Output ($BB)']*1e6)*0.23*0.27/10
}






################################################################################
################################################################################
# 3. Make Exposure plots
################################################################################
################################################################################

# Set plot data (options: Raw or Projected)
df.WIOD.plot = df.WIOD.exp.best[df.WIOD.exp.best$Year==2014 & df.WIOD.exp.best$Code!="U" & df.WIOD.exp.best$Code!="TOTII_new",]


df.WIOD.plot = df.WIOD.plot[order(df.WIOD.plot$expshare_automation_binary_laborsharegdp_weight),]
df.WIOD.plot['Code_f'] = factor(df.WIOD.plot$Code,levels=df.WIOD.plot$Code)
df.WIOD.plot['Description_short'] = df.WIOD.plot['Description']
df.WIOD.plot['Description_f2'] = 0
for(i in c(1:dim(df.WIOD.plot)[1])){
  df.WIOD.plot[i,'Description_f2'] = i
  if(nchar(df.WIOD.plot[i,"Description_short"])>43) {
    df.WIOD.plot[i,"Description_short"] = paste0(substr(df.WIOD.plot[i,"Description_short"],1,40),"...")
  }
  
}
df.WIOD.plot['Description_f'] = factor(df.WIOD.plot$Description_short,levels=df.WIOD.plot$Description_short)




pdf("TablesAndFigures/ISICRev4_Industry_Wagebillweighted_Exposure.pdf",width=3,height=5)
ggplot(df.WIOD.plot[order(df.WIOD.plot$expshare_automation_binary_laborsharegdp_weight),]) +
  geom_bar(aes(x=expshare_automation_binary_laborsharegdp_weight*100,y=Description_f,fill=Description_f2),stat='identity',width=.5) +
  theme_classic() + theme(axis.text.y = element_text(size = 5),axis.text.x = element_text(size = 8),text = element_text(size=8)) +
  scale_fill_gradient(low = "lightgreen", high = "darkgreen") +
  theme(legend.position = "none") +
  scale_x_continuous(position = "top",breaks=c(0,0.1,0.2),labels=c("0%","0.1%","0.2%")) +
  xlab(expression(frac(Delta*A[k],A[k])*"")) + ylab("ISIC Rev.4 Industry") #+ ggtitle("Industry wage-bill weighted exposure shares (ISIC Rev.4)")
dev.off()



# Plot energy intensity and emissions intensity
ggplot() +
  geom_point(data = df.WIOD.plot,aes(x=Description_f,y=`Energy Intensity (PJ/$BB)`)) +
  geom_hline(yintercept = sum(df.WIOD.plot$`Total Energy Use (PJ)`)/sum(df.WIOD.plot$`Output ($BB)`),color="grey",linetype="dashed") +
  theme_classic() +
  theme(axis.text.x=element_text(angle=90, vjust = 0, hjust=1),axis.title.x=element_blank()) + 
  scale_y_continuous(expand = c(0, 0), limits = c(0, NA))

ggplot() +
  geom_point(data = df.WIOD.plot,aes(x=Description_f,y=`Emissions Intensity (kt/PJ)`)) +
  geom_hline(yintercept = sum(df.WIOD.plot$`CO2 Emissions (kt)`)/sum(df.WIOD.plot$`Total Energy Use (PJ)`),color="grey",linetype="dashed") +
  theme_classic() +
  theme(axis.text.x=element_text(angle=90, vjust = 0, hjust=1),axis.title.x=element_blank()) + 
  scale_y_continuous(expand = c(0, 0), limits = c(0, NA))


################################################################################
################################################################################
# 4. Make Energy and Emissions change plots (Aggregate)
################################################################################
################################################################################
# Set exposure score to be used
expscore = "automation_binary"
expscore.weight = "_laborsharegdp_weight"


df.WIOD.plot['Change in GDP'] = (df.WIOD.plot[paste0('expshare_',expscore,expscore.weight)])*df.WIOD.plot$`Output ($BB)`
df.WIOD.plot['Change in Energy'] = df.WIOD.plot['Change in GDP']*df.WIOD.plot$`Energy Intensity (PJ/$BB)`
df.WIOD.plot['Change in Emissions'] = df.WIOD.plot['Change in Energy']*df.WIOD.plot['Emissions Intensity (kt/PJ)']
df.WIOD.plot['Change in GDP (upper)'] = (df.WIOD.plot[paste0('expshare_',expscore,"_upper",expscore.weight)])*df.WIOD.plot$`Output ($BB)`
df.WIOD.plot['Change in Energy (upper)'] = df.WIOD.plot['Change in GDP (upper)']*df.WIOD.plot$`Energy Intensity (PJ/$BB)`
df.WIOD.plot['Change in Emissions (upper)'] = df.WIOD.plot['Change in Energy (upper)']*df.WIOD.plot['Emissions Intensity (kt/PJ)']
df.WIOD.plot['Change in GDP (lower)'] = (df.WIOD.plot[paste0('expshare_',expscore,"_lower",expscore.weight)])*df.WIOD.plot$`Output ($BB)`
df.WIOD.plot['Change in Energy (lower)'] = df.WIOD.plot['Change in GDP (lower)']*df.WIOD.plot$`Energy Intensity (PJ/$BB)`
df.WIOD.plot['Change in Emissions (lower)'] = df.WIOD.plot['Change in Energy (lower)']*df.WIOD.plot['Emissions Intensity (kt/PJ)']


print(paste0("Labor share GDP weighted exposure (aggregate) = ",as.character(sum((df.WIOD.plot[paste0('expshare_',expscore,"_upper",expscore.weight)])))))
# print(paste0("Change in energy use (aggregate) = ",as.character(round(sum((df.WIOD.plot[paste0('expshare_',expscore,expscore.weight)]*0.23*0.27)/10)*df.WIOD[df.WIOD$Code=="TOTII_new" & df.WIOD$Year==2014,"Output ($BB)"]*df.WIOD[df.WIOD$Code=="TOTII_new" & df.WIOD$Year==2014,"Energy Intensity (PJ/$BB)"]),3)))
print(paste0("Change in energy use (by industry) = ",as.character(round(sum(df.WIOD.plot$`Change in Energy`),3))))
print(paste0("Change in emissions (by industry) = ",as.character(round(sum(df.WIOD.plot$`Change in Emissions`),3))))
print(paste0("Change in energy use (aggregate) = ",as.character(round(sum(df.WIOD.plot['Change in GDP'])*sum(df.WIOD.plot$`Total Energy Use (PJ)`)/sum(df.WIOD.plot$`Output ($BB)`),3))))
print(paste0("Change in emissions (aggregate) = ",as.character(round(sum(df.WIOD.plot['Change in GDP'])*sum(df.WIOD.plot$`Total Energy Use (PJ)`)/sum(df.WIOD.plot$`Output ($BB)`)*sum(df.WIOD.plot$`CO2 Emissions (kt)`)/sum(df.WIOD.plot$`Total Energy Use (PJ)`),3))))

# Export relevant data as table
write.csv(df.WIOD.plot,"Data/Output/FinalEstimates.csv",row.names=F)

# Set codes to highlight (Publishing activities, Education,Wholesale trade and repair of motor vehicles)
codesofinterest = c("J58","P85","G45")
# codesofinterest.color <- c("J58" = "#F46036", "P85" = "#1B998B", "G45" = "#C5D86D")
codesofinterest.color <- c("J58" = "red", "P85" = "blue", "G45" = "green")



# Create a data frame

xvar = df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,"Output ($BB)"]
yvar_min = df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,"expshare_automation_binary_lower_laborsharegdp_weight"]
yvar_max = df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,"expshare_automation_binary_upper_laborsharegdp_weight"]
x_range = seq(min(xvar),max(xvar),(max(xvar)-min(xvar))/100)
df001 = data.frame("x"=x_range,"y"=1e-5/(x_range))
df10 = data.frame("x"=x_range,"y"=1e-3/(x_range))
df100 = data.frame("x"=x_range,"y"=1e-1/(x_range))
df10000 = data.frame("x"=x_range,"y"=1e1/(x_range))
gdp = sum(df.WIOD.plot$`Output ($BB)`)
plt.ratio = 2/5
x.pos = exp(log(min(x_range)+0.00025*(max(x_range)-min(x_range))))

plt1 <- ggplot() +
  geom_point(data=df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,],aes(`Output ($BB)`,expshare_automation_binary_laborsharegdp_weight),color="grey") +
  geom_point(data=df.WIOD.plot[df.WIOD.plot$Code %in% codesofinterest,],aes(`Output ($BB)`,expshare_automation_binary_laborsharegdp_weight,color=Code)) +
  geom_linerange(data=df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,],aes(x=`Output ($BB)`,ymin=expshare_automation_binary_lower_laborsharegdp_weight,ymax=expshare_automation_binary_upper_laborsharegdp_weight),color="grey") +
  geom_linerange(data=df.WIOD.plot[df.WIOD.plot$Code %in% codesofinterest,],aes(x=`Output ($BB)`,ymin=expshare_automation_binary_lower_laborsharegdp_weight,ymax=expshare_automation_binary_upper_laborsharegdp_weight,color=Code)) +
  scale_x_continuous(trans='log10',breaks=c(30,100,300,1000,3000),labels=c("30","100","300","1,000","3,000")) + scale_y_continuous(trans='log10',limits=c(1e-7,2e-2)) +
  geom_line(data = df001,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=x.pos,y=1e-6,label="$10 Thousand",angle=180/pi*atan(-1*plt.ratio),hjust = 0) +
  geom_line(data = df10,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=x.pos,y=1e-4,label="$1 Million",angle=180/pi*atan(-1*plt.ratio),hjust = 0) +
  geom_line(data = df100,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=x.pos,y=1e-2,label="$100 Million",angle=180/pi*atan(-1*plt.ratio),hjust = 0) +
  geom_line(data = df10000,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=750,y=2e-2,label="$10 Billion",angle=180/pi*atan(-1*plt.ratio),hjust = 0) +
  scale_color_manual(values = codesofinterest.color) +
  ylab(expression("Change in productivity (" * Delta * A[k] / A[k] * ", %)")) + xlab(expression("Economic output (" * y[k] * ", $B)")) +
  theme_classic(base_size=15) + theme(legend.position="none") #+ coord_fixed(ratio=plt.ratio)

print(plt1)


pdf("TablesAndFigures/ShockvsValueAdded.pdf",width=5,height=5)
print(plt1)
dev.off()


xvar = df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,"Energy Intensity (PJ/$BB)"]
yvar_min = df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,"Change in GDP (lower)"]
yvar_max = df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,"Change in GDP (upper)"]
x_range = seq(min(xvar),max(xvar),(max(xvar)-min(xvar))/1000)
df0001 = data.frame("x"=x_range,"y"=1e-5/x_range)
df10 = data.frame("x"=x_range,"y"=1e-2/x_range)
df100 = data.frame("x"=x_range,"y"=1/x_range)
df10000 = data.frame("x"=x_range,"y"=1e2/x_range)
plt.ratio = 2/3
x.pos = exp(log(min(x_range)+0.0000025*(max(x_range)-min(x_range))))

plt2 <- ggplot() +
  geom_point(data=df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,],aes((`Energy Intensity (PJ/$BB)`),`Change in GDP`),color="grey") +
  geom_point(data=df.WIOD.plot[df.WIOD.plot$Code %in% codesofinterest,],aes((`Energy Intensity (PJ/$BB)`),`Change in GDP`,color=Code)) +
  geom_linerange(data=df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,],aes(x=(`Energy Intensity (PJ/$BB)`),ymin=`Change in GDP (lower)`,ymax=`Change in GDP (upper)`),color="grey") +
  geom_linerange(data=df.WIOD.plot[df.WIOD.plot$Code %in% codesofinterest,],aes(x=(`Energy Intensity (PJ/$BB)`),ymin=`Change in GDP (lower)`,ymax=`Change in GDP (upper)`,color=Code)) +
  geom_line(data = df0001,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=x.pos,y=3e-3,label="10 GJ",angle=180/pi*atan(-1*plt.ratio),hjust=0) +
  geom_line(data = df10,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=x.pos,y=3,label="10 TJ",angle=180/pi*atan(-1*plt.ratio),hjust=0) +
  geom_line(data = df100,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=1.1e-1,y=1.5e1,label="1 PJ",angle=180/pi*atan(-1*plt.ratio),hjust=0) +
  geom_line(data = df10000,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=0.8e1,y=2e1,label="100 PJ",angle=180/pi*atan(-1*plt.ratio),hjust=0) +
  scale_color_manual(values = codesofinterest.color) +
  scale_x_continuous(trans='log10',breaks=c(1e-2,1e-1,1,1e1,1e2),labels=c("0.01","0.1","1","10","100")) + scale_y_continuous(trans='log10',breaks=c(1e-5,1e-3,1e-1,1e1),labels=c("0.01","1","100","10,000"),limits=c(0.5e-5,2e1)) +#,breaks=c(1e-3,1,1e3),limits = c(1e-3,1e6)) +
  ylab(expression("Change in economic output (" * Delta * y[k] * ", $M)")) + xlab(expression("Energy intensity (" * nu[k] * ", TJ/$M)")) +
  theme_classic(base_size=15) + theme(legend.position="none") #+ coord_fixed(ratio=plt.ratio)

print(plt2)

pdf("TablesAndFigures/ChangeinGDPvsEnergyIntensity.pdf",width=5,height=5)
print(plt2)
dev.off()

xvar = df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,"Emissions Intensity (kt/PJ)"]
yvar_min = df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,"Change in Energy (lower)"]
yvar_max = df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,"Change in Energy (upper)"]
x_range = seq(min(xvar),max(xvar),(max(xvar)-min(xvar))/100)
y_range = 1/x_range
df00001 = data.frame("x"=x_range,"y"=1e-3/x_range)
df1 = data.frame("x"=x_range,"y"=1e-1/x_range)
df100 = data.frame("x"=x_range,"y"=1e1/x_range)
df1000 = data.frame("x"=x_range,"y"=1e3/x_range)
plt.ratio = 2/8
x.pos = exp(log(min(x_range)+0.0000025*(max(x_range)-min(x_range))))

plt3 <- ggplot() +
  geom_point(data=df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,],aes(`Emissions Intensity (kt/PJ)`,`Change in Energy`),color="grey") +
  geom_point(data=df.WIOD.plot[df.WIOD.plot$Code %in% codesofinterest,],aes(`Emissions Intensity (kt/PJ)`,`Change in Energy`,color=Code)) +
  geom_linerange(data=df.WIOD.plot[df.WIOD.plot['expshare_automation_binary_laborsharegdp_weight']!=0,],aes(x=(`Emissions Intensity (kt/PJ)`),ymin=`Change in Energy (lower)`,ymax=`Change in Energy (upper)`),color="grey") +
  geom_linerange(data=df.WIOD.plot[df.WIOD.plot$Code %in% codesofinterest,],aes(x=(`Emissions Intensity (kt/PJ)`),ymin=`Change in Energy (lower)`,ymax=`Change in Energy (upper)`,color=Code)) +
  scale_x_continuous(trans='log10',breaks=c(3,10,30,100),labels=c(3,10,30,100)) + scale_y_continuous(trans='log10',limits = c(0.5e-5,2e2),breaks=c(1e-4,1e-2,1,1e2),labels=c("0.1","10","1,000","100,000")) + #,breaks=c(1e-3,1,1e3),limits = c(1e-5,1e5)) +
  geom_line(data = df00001,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=x.pos,y=0.5e-3,label=expression("1 tC"*O[2]),angle=180/pi*atan(-1*plt.ratio),hjust=0) +
  geom_line(data = df1,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=x.pos,y=0.5e-1,label=expression("100 tC"*O[2]),angle=180/pi*atan(-1*plt.ratio),hjust=0) +
  geom_line(data = df100,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=x.pos,y=0.5e1,label=expression("10 ktC"*O[2]),angle=180/pi*atan(-1*plt.ratio),hjust=0) +
  geom_line(data = df1000,aes(x=x,y=y),linewidth=0.5,color="grey",linetype="dashed") +
  annotate(geom="text",x=9,y=2e2,label=expression("1 MtC"*O[2]),angle=180/pi*atan(-1*plt.ratio),hjust=0) +
  scale_color_manual(values = codesofinterest.color) +
  ylab(expression("Change in energy (" * Delta * E[k] * ", TJ)")) + xlab(expression("Emissions intensity (" * mu[k] * ", " * tCO[2] * "/TJ)")) +
  theme_classic(base_size=15) + theme(legend.position="none") #+ coord_fixed(ratio=plt.ratio)

print(plt3)


pdf("TablesAndFigures/ChangeinEnergyvsEmissionsIntensity.pdf",width=5,height=5)
print(plt3)
dev.off()


grid.arrange(plt1, plt2, plt3, ncol=3)





################################################################################
################################################################################
# 5. Make Energy and Emissions change plots (By Industry)
################################################################################
################################################################################
df.WIOD.plot <- df.WIOD.plot[order(df.WIOD.plot$Code),]
df.WIOD.plot['Description_f'] = factor(df.WIOD.plot$Description_short,levels=df.WIOD.plot$Description_short)

# Plot impact on GDP by Industry
pdf("TablesAndFigures/Industry_GDP_Impact.pdf",width=3,height=5)
ggplot() +
  geom_bar(data=df.WIOD.plot[order(df.WIOD.plot$Code),],aes(x=`Change in GDP`,y=Description_f,fill=Description_f2),stat='identity',width=.5) +
  theme_classic() + theme(axis.text.y = element_text(size = 5),axis.text.x = element_text(size = 5),text = element_text(size=8)) +
  scale_fill_gradient(low = "lightgreen", high = "darkgreen") +
  theme(legend.position = "none") +
  scale_x_continuous(position = "top") +
  xlab("Change in GDP ($BB)") + ylab("ISIC Rev.4 Industry") 
dev.off()


# Plot impact on Energy Use by Industry
pdf("TablesAndFigures/Industry_EnergyUse_Impact.pdf",width=3,height=5)
ggplot() +
  geom_bar(data=df.WIOD.plot,aes(x=`Change in Energy`,y=Description_f,fill=Description_f2),stat='identity',width=.5) +
  theme_classic() + theme(axis.text.y = element_text(size = 5),axis.text.x = element_text(size = 5),text = element_text(size=8)) +
  scale_fill_gradient(low = "lightgreen", high = "darkgreen") +
  theme(legend.position = "none") +
  scale_x_continuous(position = "top") +
  xlab("Change in Energy (PJ)") + ylab("ISIC Rev.4 Industry") #+ ggtitle("Industry wage-bill weighted exposure shares (ISIC Rev.4)")
dev.off()


# Plot impact on Emissions by Industry
pdf("TablesAndFigures/Industry_Emissions_Impact.pdf",width=3,height=5)
ggplot() +
  geom_bar(data=df.WIOD.plot,aes(x=`Change in Emissions`,y=Description_f,fill=Description_f2),stat='identity',width=.5) +
  theme_classic() + theme(axis.text.y = element_text(size = 5),axis.text.x = element_text(size = 5),text = element_text(size=8)) +
  scale_fill_gradient(low = "lightgreen", high = "darkgreen") +
  theme(legend.position = "none") +
  scale_x_continuous(position = "top") +
  xlab("Change in Emissions (ktCO2)") + ylab("ISIC Rev.4 Industry") #+ ggtitle("Industry wage-bill weighted exposure shares (ISIC Rev.4)")
dev.off()



################################################################################
################################################################################
# 6. Make sensitivity plots
################################################################################
################################################################################


x_range = seq(0.00001,1,0.05)
y_range = seq(0.00001,1,0.05)
grid = expand.grid("x"=x_range,"y"=y_range)
grid$z = grid$x*grid$y*sum(df.WIOD.plot$`Change in Energy`/0.27/0.23)*10
Energy.change = sum(df.WIOD.plot$`Change in Energy`)

df.agg = data.frame("cost savings"=0.27,"profitability"=0.23)

sensplt1 <- ggplot() +
  geom_contour_filled(data = grid,aes(x=x,y=y,z=z,fill=..level..),breaks=c(0,100,Energy.change*10,1000,3000,1000000)) +
  geom_contour(data = grid,aes(x=x,y=y,z=z),color="black",breaks=c(100,Energy.change*10,1000,3000)) +
  geom_point(data=df.agg,aes(cost.savings,profitability),size=3) +
  annotate(geom="text",x=0.065,y=0.5,label="10 PJ",angle=-80) +
  annotate(geom="text",x=0.135,y=0.55,label="28 PJ",angle=-76) +
  annotate(geom="text",x=0.35,y=0.68,label="100 PJ",angle=-61) +
  annotate(geom="text",x=0.795,y=0.87,label="300 PJ",angle=-44) +
  scale_fill_brewer(palette = "Blues") +
  scale_x_continuous(limits=c(0,1),breaks=c(0,0.25,0.50,0.75,1),labels=c("0%","25%","50%","75%","100%"),expand = c(0, 0)) + scale_y_continuous(limits=c(0,1),breaks=c(0,0.25,0.50,0.75,1),labels=c("0%","25%","50%","75%","100%"),expand = c(0, 0)) +
  ylab("Profitably exposed tasks") + xlab("Avg. labor cost savings") +
  theme_classic(base_size=15) + theme(legend.position="none") + expand_limits(x = 0, y = 0) + theme(plot.margin=margin(t=20,r=20)) + guides(fill=guide_legend(title=expression(Delta*"Energy Use (PJ)")))

print(sensplt1)



x_range = seq(0.00001,1,0.05)
y_range = seq(0.00001,1,0.05)
grid = expand.grid("x"=x_range,"y"=y_range)
grid$z = grid$x*grid$y*sum(df.WIOD.plot$`Change in Emissions`/0.27/0.23)*10
emissions.change = sum(df.WIOD.plot$`Change in Emissions`)*10

df.agg = data.frame("cost savings"=0.27,"profitability"=0.23)

sensplt2 <- ggplot() +
  geom_contour_filled(data = grid,aes(x=x,y=y,z=z,fill=..level..),breaks=c(0,3000,emissions.change,30000,60000,1000000)) +
  geom_contour(data = grid,aes(x=x,y=y,z=z),color="black",breaks=c(3000,emissions.change,30000,60000)) +
  geom_point(data=df.agg,aes(cost.savings,profitability),size=3) + 
  annotate(geom="text",x=0.065,y=0.5,label=expression("300 ktC"*O[2]),angle=-81) +
  annotate(geom="text",x=0.14,y=0.55,label=expression("900 ktC"*O[2]),angle=-73) +
  annotate(geom="text",x=0.345,y=0.66,label=expression("3 MtC"*O[2]),angle=-60) +
  annotate(geom="text",x=0.565,y=0.78,label=expression("6 MtC"*O[2]),angle=-50) +
  scale_fill_brewer(palette = "Oranges") +
  scale_x_continuous(limits=c(0,1),breaks=c(0,0.25,0.50,0.75,1),labels=c("0%","25%","50%","75%","100%"),expand = c(0, 0)) + scale_y_continuous(limits=c(0,1),breaks=c(0,0.25,0.50,0.75,1),labels=c("0%","25%","50%","75%","100%"),expand = c(0, 0)) +
  ylab("Profitably exposed tasks") + xlab("Avg. labor cost savings") +
  theme_classic(base_size=15) + theme(legend.position="none") + expand_limits(x = 0, y = 0) + theme(plot.margin=margin(t=20,r=20)) + guides(fill=guide_legend(title=expression(Delta*"Emissions (ktC"*O[2]*")")))

print(sensplt2)

grid.arrange(sensplt1,sensplt2,ncol=2)


pdf("TablesAndFigures/ChangeinEnergySensitivity.pdf",width=5,height=5)
print(sensplt1)
dev.off()


pdf("TablesAndFigures/ChangeinEmissionsSensitivity.pdf",width=5,height=5)
print(sensplt2)
dev.off()






